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PREFACE 


This  thesis  compares  the  results  of  using  the  trigonom- 
etric and  conventional  approaches  to  the  finite  difference 
calculus  incorporating  both  virtual  work  and  equilibrium 
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ABSTRACT 

A relatively  new  trigonometric  approach  to  the  finite 
difference  calculus  was  applied  to  the  problem  of  beam  buckling 
as  represented  by  both  virtual  work  and  equilibrium  equations. 
The  trigonometric  functions  were  varied  by  adjusting  a wave- 
length parameter  in  the  approximating  Fourier  series.  Values 
of  the  critical  force  obtained  from  the  modified  approach  for 
beams  with  a variety  of  boundary  conditions  were  compared  to 
results  using  the  conventional  finite  difference  method.  The 
trigonometric  approach  produced  significantly  more  accurate 
approximations  for  the  critical  force  than  the  conventional 
approach  for  a relatively  wide  range  in  values  of  the  wave- 
length parameter;  and  the  optimizing  value  of  the  wavelength 
parameter  corresponded  to  the  half  wavelength  of  the  buckled 
mode  shape.  It  was  found  from  a modal  analysis  that  the  most 
accurate  solutions  are  obtained  when  the  approximating  function 
closely  represents  the  actual  displacement  function.  It  is 
more  difficult  to  select  a satisfactory  value  of  the j^velength 
parameter  for  the  equilibrium  equation  which  makes  the  Virtual 
work  equation  more  attractive  for  practical  applications.  The 
buckled  mode  shape  (or  eigenfunction)  is  predicted  with  high 
accuracy  regardless  of  the  value  of  the  wavelength  parameter. 

A comparison  of  the  virtual  work  and  the  Galerkin  approaches 
identified  marked  similarities  between  the  two  methods. 
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APPLICATION  OF  TRIGONOMETRIC  AND  CONVENTIONAL 
FINITE  DIFFERENCE  APPROXIMATIONS 
TO  BEAM  BUCKLING 

I.  Introduction 


Background 

One  of  the  greatest  challenges  facing  the  space  age  engin- 
eer involves  solving  complex  differential  equations  with  a 
high  degree  of  accuracy  and  efficiency.  An  exact  solution  to 
most  modern  day  problems  is  possible  only  if  major  simplifying 
assumptions  are  made  which  linearize  the  differential  equations. 
These  analytical  solutions  are  useful  in  providing  guidelines 
for  the  engineer;  but  due  to  the  simplifications  involved, 
they  do  not  provide  the  degree  of  accuracy  required  for  the 
majority  of  engineering  applications.  The  advent  of  the  high 
speed  computer  made  it  possible,  for  the  first  time,  to  obtain 
solutions  to  complex  mathematical  models  which  realistically 
represented  the  actual,  non-linear  world.  The  accuracy  of 
these  numerical  solutions  is  limited  only  by  the  precision  of 
the  numerical  technique  and  by  the  designer's  ability  to  con- 
struct a realistic  mathematical  model.  As  a result,  primary 
emphasis  in  recent  years  has  been  placed  on  the  development  of 
new  numerical  schemes  which  can  solve  more  sophisticated  prob- 
lems at  a faster  rate  and  with  increased  accuracy. 

The  finite  difference  calculus  is  a numerical  technique 
for  solving  differential  equations  which  (like  most  other 
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numerical  techniques)  divides  the  continuum  into  a finite 
number  of  degrees  of  freedom.  The  conventional  technique 
derives  algebraic  approximations  for  the  derivatives  of  a 
function  by  combining  various  forms  of  the  Taylor  series  ex- 
pansion. A relatively  new  technique  developed  by  M.  Stein 
and  J.  Housner  [1]  uses  trigonometric  expressions  to  approxi- 
mate the  derivatives  in  a differential  equation.  The  trigon- 
ometric expressions  are  derived  from  the  truncated  Fourier 
series  instead  of  the  Taylor  series  and  are  believed  to  provide 
a closer  approximation  for  functions  with  sinusoidal  charac- 
teristics. The  succeeding  sections  of  this  thesis  compare 
the  conventional  and  trigonometric  approaches  to  the  finite 
difference  calculus  when  applied  to  the  problem  of  beam  buck- 
ling. 

Beam  buckling  is  one  of  the  most  severe  problems  in 
structural  design,  and  the  prediction  of  the  critical,  or 
buckling,  force  is  of  primary  concern.  For  example,  the 
landing  gear  strut  of  an  aircraft  during  landing  will  be  sub- 
jected to  relatively  large  axial  forces.  The  minimum  force 
which  will  cause  buckling  is  a primary  determinant  of  strut 
design.  When  this  critical  axial  force  is  applied  to  a beam, 
sudden  bending  occurs  because  the  beam  is  in  a form  of  unstable 
equilibrium.  Therefore,  the  critical  load  can  be  defined  as 
that  axial  force  which  will  just  maintain  a slightly  deflected 
configuration  [2].  Forces  above  the  critical  load  will  result 
in  an  acceleration  of  the  rate  of  deflection. 

Numerous  methods  have  been  used  to  calculate  the  critical 
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force.  Analytical  solutions  can  be  obtained  for  simplified 
problems  involving  beams  with  uniform  cross  section  and  sim- 
ple boundary  conditions  [3-7].  More  complex  problems  have 
been  solved  by  a variety  of  techniques  including  Ral ei gh-Ri tz , | 

Galerkin,  finite  element,  and  use  of  Fourier  series  [3,8-11]. 

These  techniques  are  useful  but  tend  to  be  lengthy  and  tedious.  | 

Application  of  the  finite  difference  method  has  also  been  dis- 
cussed at  length  in  the  literature  [12-16].  This  approach 
provides  the  solution  for  beams  with  complex  boundary  condi- 
tions and  variable  cross  section  as  easily  as  uniform,  simple 
beams.  Most  of  the  previous  research  in  the  area  of  finite 
differences  deals  exclusively  with  beams  which  have  either 
simple  or  clamped  boundary  conditions. 

Purpose 

The  purpose  of  this  thesis  is  to  compare  the  trigonome- 
tric and  conventional  approaches  to  the  finite  difference 
method  in  calculating  the  critical  force  for  a one  dimensional 
beam.  Solutions  were  obtained  for  both  the  equilibrium  and 
virtual  work  equations  using  a functional  parameter  based  upon 
the  structure's  buckle  mode.  The  criteria  for  comparison  are 
accuracy  and  efficiency. 

General  Approach 

Since  Stein  and  Housner  dealt  exclusively  with  the  second 
order  virtual  work  equation  during  their  investigation  of  the 
trigonometric  approach  to  the  finite  difference  method,  they 
made  no  attempt  to  develop  trigonometric  approximations  for 


3 


third  and  higher  order  derivatives  [1].  Consequently,  the 
first  step  of  this  research  was  to  derive  the  trigonometric 
finite  difference  expression  for  the  fourth  derivative.  This 
derivation  provided  the  remaining  tool  needed  to  evaluate 
the  fourth  order  equilibrium  equation  using  both  the  trigon- 
ometric and  conventional  (polynomial)  approaches.  To  accom- 
plish this  task,  the  domain  of  the  beam  was  separated  into  a 
set  of  node  points.  Finite  difference  approximations  were 
developed  at  each  node  point,  and  the  results  for  each  node 
were  substituted  Into  the  equilibrium  equation.  This  substitu- 
tion resulted  In  an  eigenvalue  problem  from  which  the  critical 
force  was  solved.  Various  combinations  of  four  end  restraints 
were  considered.  These  restraints  consisted  of  pinned, 
clamped,  guided,  and  free. 

A similar  procedure  was  conducted  for  the  virtual  work 
approach.  The  finite  difference  expressions  were  substituted 
into  the  virtual  work  equation  and  integration  was  performed 
over  the  length  of  the  beam  using  the  trapezoid  rule.  A sort- 
ing technique  was  used  to  isolate  the  virtual  displacements 
and  establish  a set  of  algebraic  equations  leading  to  an 
eigenvalue  problem.  The  critical  force  was  obtained  from  the 
solution  of  this  eigenvalue  problem. 

Previous  work  concentrated  on  applying  the  trigonometric 
approach  to  the  buckling  of  plates  with  simple  and  clamped 
boundary  conditions  [1,16].  Solutions  were  obtained  for  the 
virtual  work  equation  only.  This  thesis  deals  with  the  problem 
of  beam  buckling  and  investigates  a wider  range  of  boundary 
conditions.  Solutions  to  both  the  equilibrium  and  virtual  work 
equations  were  obtained.  4 


1 1 . Theory 


Assumptions 

The  equations  for  beam  buckl  1 ng  were  derived  with  the 
following  assumptions: 

1.  Sections  of  the  beam  normal  to  the  longitudinal  axis 
remain  plane  during  buckling. 

2.  The  cross  sectional  area  of  the  beam  Is  small  com- 
pared to  the  beam  length. 

3.  Displacements  remain  Infinitesimally  small. 

4.  The  beam  is  perfectly  flat  prior  to  buckling. 

5.  The  beam  Is  composed  of  a homogeneous.  Isotropic 
material . 

6.  The  relationship  for  strain  energy  used  to  derive  the 
virtual  work  equation  can  be  approximated  by 

“ ■ lllv  ? “x'* 

where  e is  the  strain  in  the  axial  direction  and  a Is  the 

X X 

stress  In  the  axial  direction  [5]. 

Equi 1 1 bri urn  Pi ff erentlal  Equati on 

The  development  of  equations  of  equilibrium  In  elastic 
theory  date  back  to  Poisson  and  Cauchy  in  the  nineteenth  cen- 
tury [6].  The  equilibrium  differential  equation  describing 
the  deflection  of  a beam  subjected  to  a constant  axial  load 
may  be  expressed  as 


where  V is  the  out  of  plane  displacement  and  the  modulus  of 
elasticity,  E,  and  the  moment  of  inertia,  I,  are  assumed  to  be 
constant.  The  coordinate  system  and  placement  of  the  axial 
force,  P,  are  shown  in  Fig.  2.1. 


The  derivation  of  Eq  (2-1)  can  be  found  in  numerous  texts  and 
is  presented  in  Appendix  A [7,5,17]. 

Virtual  Work  Equation 

The  virtual  work  principle  can  be  stated  as  follows.  If 
a system  acted  upon  by  a number  of  internal  and  external  for- 
ces is  in  equilibrium,  then  the  total  work  done  by  these 
forces  due  to  a small  virtual  displacement  is  zero  [47].  This 
axiom  can  be  restated  as 

6W  = -6W.  (2-3) 

e T 

where  6W  and  6W.  represent  the  external  and  internal  work  due 
e 1 

to  a small  virtual  displacement.  Internal  work  is  related  to 
potential  energy,  U,  by  W^  = -U.  The  virtual  work  equation 
for  a one  dimensional  beam  with  an  axial  force  is 
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where  V represents  actual  displacements  and  6V  represents  vir- 
tual displacements  [19].  The  derivation  of  the  virtual  work 
equation  is  shown  In  Appendix  B. 

Boundary  Conditions 

The  following  end  constraints  were  applied  to  the  beam  In 
various  combinations:  pinned,  clamped,  guided,  and  free.  These 
constraints  were  specifically  chosen  so  that  the  full  range  of 
boundary  conditions  could  be  explored.  The  following  relation- 
ships exist  between  the  end  constraints  and  the  mathematical 
representation  of  the  boundary  condition: 


1. 

Pinned : 

V(0)  = 0 

V’'(0) 

= 0 

2. 

Cl amped : 

V(0)  = 0 

V'(0) 

= 0 

3. 

Gul ded : 

V'(0)  = 0 

V"  (0) 

+ Tv  (0)  = 0 

4. 

Free : 

V"(0)  = 0 

V"  (0) 

+ Tv(o)  = 0 

where  P is  a non-dimensional  load  parameter  defined  by  the 
following  expression:  P = PL  /El 
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III.  Numerical  Technique 


Finite  Difference  Cal cul us 

Conventional  Approach.  The  finite  difference  calculus  is 
a numerical  technique  for  solving  differential  equations  which 
divides  a continuous  function  into  a finite  number  of  degrees 
of  freedom.  This  method  was  introduced  by  Boole  and  others  in 
the  nineteenth  century  [20].  Since  that  time  finite  differen- 
ces have  been  given  a great  deal  of  attention  as  a vehicle  for 
solving  non-linear  differential  equations  in  a variety  of  eng- 
ineering disciplines  [21-24].  As  mentioned  in  the  introduction 
to  this  report,  the  finite  difference  technique  i s particul arl y 
well-suited  to  the  field  of  structural  engineering  in  solving 
buckling  problems. 

The  concept  underlying  the  finite  difference  calculus  is 
relatively  simple.  The  function  under  investigation  is  first 
separated  into  a number  of  grid  points.  Fig.  3.1  represents 
a function,  v = f(x),  which  for  this  analysis  can  be  considered 
to  be  the  unknown  deflections  of  a beam. 
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Fig.  3,1.  Finite  Difference  Grid 
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Grid  points  (indicated  by  circles  in  Fig.  3.1)  represent  dis- 
crete points  at  which  the  function  is  defined.  These  points 
are  separated  by  the  distance  h.  Reference  points  (Indicated 
by  an  X)  refer  to  points  at  which  the  finite  difference  approx- 
imations are  desired.  Grid  points  and  reference  points  may 
coincide,  but  this  Is  not  mandatory. 

The  conventional  approach  to  the  finite  difference  calcu- 
lus derives  algebraic  approximations  for  the  derivatives  In  a 
differential  equation  by  combining  various  forms  of  the  trunca- 
ted Taylor  series.  The  basic  form  of  the  Taylor  series  Is 

V(X+h)  = V(X)  + hV'(X)  + |h^  V"(X)  + ^h^  V"  (X) 

+ ^ V^''(X)  + •••  (3-1) 

In  short  hand  notation  this  expression  can  be  written  as 

''vl  “ ^0  ^ ^ ''o  ^ (3-2) 

where  subscripts  are  used  to  indicate  the  displacement  from 
the  reference  point.  Similarly,  the  Taylor  series  expansion 
for  a small  displacement  to  the  left  of  the  reference  point  Is 


'-1 


= V 


hV 


ih2 


- 


6^ 


Y III 
0 


J 

24 


4 „1v 


+ V 


(3-3) 


By  subtracting  Eq  (3-2)  from  Eq  (3-3)  and  rearranging  terms, 
the  following  expression  Is  obtained 

’i  ■ iV  - (’r’ ’o  *•••  (’-‘'l 

2 

The  term  (h  represents  the  error  approximation  and  the 
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remainder  of  Eq  (3-4)  Is  the  conventional  finite  difference 
approximation  for  the  first  derivative: 


(3-5) 

If  Eq  (3-2)  Is  added  to  Eq  (3-3),  the  finite  difference  ex- 
pression for  the  second  derivative  Is  obtained: 

(''+l-3''o*''.l)  - T7  (3-6) 

h 

where  (h^/12)V^'^  Is  the  error  approximation. 

Conventional  finite  difference  approximations  for  higher 
order  derivatives  can  be  derived  In  a similar  manner.  The 
development  of  the  fourth  derivative  expression  requires  one 
additional  Taylor  series  expansion  In  each  direction  resulting 
In 

''i’  ■ jT  (''.2-^''-1*6Vo-4V^,*V,2)  (3-7) 

A geometric  Interpretation  for  this  conventional  approach 
can  be  .obtained  by  passing  an  nth  degree  polynomial  through 
the  set  of  n+1  data  points  which  define  the  function  under 
Investigation.  In  fact.  It  Is  possible  to  derive  finite  diff- 
erence approximations  by  fitting  a polynomial  to  the  data 
points  [22].  For  example,  the  expression  for  the  second  deriv- 
ative can  be  found  by  fitting  the  parabola 

V(X)  = AX^  BX  + C (3-8) 

to  the  points  x = -h,0,h.  Note  that  these  points  may  be 
chosen  without  loss  of  generality.  The  second  derivative  of 
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T 


Eq  (3-8)  1s 


V"(x)  = 2A 


(3-9) 


Fitting  the  parabola  to  the  three  points  gives 

V 1 = Ah^  - Bh  + C (3-10) 

Vfl  * C (3-11) 

= Ah^  + Bh  + C (3-12) 

The  following  expression  Is  obtained  by  adding  Eqs  (3-10)  and 
(3-12): 

V_^  + = 2Ah^  + 2C  (3-13) 

Substitutions  from  Eqs  (3-9)  and  (3-11)  provide  the  equation 

V-i  + = v]^(x)h^  + 2Vj,  (3-M) 

and  solving  for  Vp(x)  yields  an  expression  for  the  second 
derivative  which  Is  Identical  to  Eq  (3-6) 

v;  = i (V.,  - 2V  * (3-15) 

n 

Trigonometric  Approach.  If  the  function  being  approxi- 
mated Is  In  fact  a polynomial,  the  conventional  finite  differ- 
ence approach  will  provide  good  approximations  for  the 
derivatives  of  the  function.  However,  a function  with  sinu- 
soidal characteristics,  such  as  the  buckling  mode  shape,  may 
be  better  approximated  by  passing  a series  of  trigonometric 
'curves  through  the  set  of  data  points  defining  the  function. 
With  this  concept  In  mind,  Stein  and  Honsner  developed  a trig- 
onometric approach  to  the  finite  difference  calculus  [1].  The 
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derivation  of  trigonometric  finite  difference  approximations 
is  similar  to  the  polynomial  derivations  with  the  exception 
that  a truncated  Fourier  series  is  used  instead  of  the  Taylor 
series. 

As  an  example  of  this  approach,  the  trigonometric  approxi- 
mation for  the  first  derivative  will  be  derived  from  the 
following  form  of  the  Fourier  series; 


V(X)  = Sin  Cos 


3r(x-x  ) tt(x-x.) 

c ..  T _g_ 


where  X represents  a variable  wavelength  parameter  [1]  and  x^^ 
is  the  reference  point.  The  derivative  of  Eq  (3-16)  with  res- 
pect to  X is 


V'(X)  . Tj  1 


Cos 


it(x-x  ) 


IT  ^(x-X  ) 

^2  I 


(3-17) 


Evaluate  Fq  (3-17)  at  x = x^  to  obtain 


V'(x„)  = f 


(3-18) 


and 


= V(x,)  J- 


(3-19) 


Now  evaluate  Eq  (3-16)  at  x^  ± h 


iT(x„+h-x„) 


iT(x„+h-x„) 


''+1  " '•'l  ■'  ’’’2  — ~ "^2  “7 — (3-20) 


Vl  •■=  Sin  f * Tj  C05  f 


-1 


T,  + T Sin  (^1)  + T Cos  (-f^-) 


(3-21) 

(3-22) 


Since  cos(-0)  = cos(C)  and  sin(-9)  - -sin(9),  Eq  (3-22)  reduces 
to 
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(3-23) 


V 


-1 


- sin  ^ + Tj  cos 


irh 

X 


Subtract  Eq  (3-23)  from  Eq  (3-21) 

- V_^  = 21^  sin  ~ (3-24) 

If  one  substitutes  Eq  (3-19)  and  rearranges  terms,  the  trigon- 
ometric approximation  for  the  first  derivative  is  obtained: 


V* 

0 


2Xsin 


"Irh 


r (V 


+ 1 


v-,) 


(3-2<) 


or 


K ‘ ; (Vi  - 

where  h = (2a/7[)  sin  (-rrh/X).  It  should  be  observed  that  Eq 
(3-26)  approaches  the  polynomial  expression  for  the  first 
derivative,  Eq  (3-5),  as  X approaches  infinity.  Their  rela- 
tionship can  be  seen  by  taking  the  limit  of  Eq  (3-25): 


lim  V'  lim  ^ 

X->-oo  ° X-v"  2Xsin 


TF  ('^+r'^-i 


(3-27) 


Since  sin(0) 
written  as 


= 0 for  small  values  of  0,  Eq  (3-27)  can  be  re- 


= (3-33) 

or 

A->’00 

Equation  (3-29)  is  identical  to  Eq  (3-5). 

Trigonometric  approximations  for  higher  order  derivatives 
can  be  derived  in  a manner  similar  to  the  derivation  of  the 
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first  derivative.  The  development  of  trigonometric  expressions 
for  the  second  and  fourth  derivatives  Is  presented  In  Appendix 
C,  and  the  results  are  Included  here  for  ease  of  reference: 


(f)'' 


1 

L (2COS9-2) 


* 16Tj  (1  - 


) 

2COS0-2)''  , 


+6V^-4V,+V2) 


(3-30) 


(3-31) 


where  Tg  Is  found  from  Eq  (C22)  In  Appendix  C and  0 = irh/X. 

(Note:  The  expression  for  was  also  derived  with  results 
similar  to  the  fourth  derivative).  It  can  again  be  shown  that 
Eqs  (3-30)  and  (3-31)  are  Identical  to  the  corresponding  poly- 
nomial expression  as  X approaches  infinity.  Proof  of  this 
characteristic  for  the  second  derivative  is  obvious.  However, 
the  fourth  derivative  is  more  difficult  to  analyze.  Each  term 
In  Eq  (3-31)  is  Indeterminate.  After  several  applications  of 
L'Hospital's  rule,  the  first  term  approaches  the  conventional 
finite  difference  approximation  for  the  fourth  derivative  and 
the  second  term  approaches  zero.  The  importance  of  this  limit- 
ing relationship  between  polynomial  and  trigonometric  expressions 
lies  in  the  fact  that  it  provides  an  indication  of  the  relative 
magnitude  of  the  error  term  for  trigonometric  approximations 
since  a great  deal  Is  known  about  the  truncation  error  for 
■conventional  finite  difference  expressions.  This  point  will 
be  discussed  further  in  succeeding  sections. 

As  indicated  previously,  the  variable  input  parameter,  X, 
is  used  to  adjust  the  wavelength  of  the  trigonometric  approx- 
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imatlng  functions.  The  determination  of  an  appropriate  value 
for  X in  Eqs  (3-25),  (3-30),  and  (3-31)  is  an  integral  part  of 
the  problem  solving  process.  Consequently,  a considerable 
amount  of  attention  is  given  in  the  remainder  of  this  thesis 
to  providing  a guide  for  choosing  the  optimizing  value  of  X. 

Half  Station  Approximations.  The  preceding  derivations 
of  finite  difference  expressions  provided  approximations  for 
derivatives  which  were  evaluated  at  the  grid  points.  It  is 
also  possible  to  derive  finite  difference  expressions  for  deriv 
atives  evaluated  midway  between  grid  points.  That  is,  the  ref- 
erence point  is  located  midway  between  two  grid  points.  This 
arrangement  is  shown  in  Fig.  3.2. 


Fig.  3.2  Half  Station  Grid  Arrangement 


The  development  of  half  station  approximations  follows  directly 
from  the  full  station  derivations.  The  Taylor  series  can  be  ex- 
pressed as  follows: 


2 

^ +i  (5-)  V"  - J (5-)^  v:  (3-32) 


2 '2'  'o  6 '0 

(I)'  vj'  .... 
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(3-33) 


'+1/2 


1 /h. 


* jV  * ■ 


Eq  (3-32)  can  be  subtracted  from  Eq  (3-33)  to  yield 


+1/2  -1/2 


-V  ,,J  - 


(3-34) 


where  h v'"  /24  represents  the  error  term.  The  error  term  for 

2 

the  full  station  approximation  was  found  to  be  h V'"  /6  in  Eq 
(3-4).  It  is  clear  that  the  first  derivative  is  more  accur- 
ately approximated  midway  between  grid  points.  On  the  other 
hand,  the  half  station  approximation  for  the  second  derivative 
is 

^ ^ v'\...(3-35) 

Comparison  of  the  error  term  in  Eq  (3-35)  with  the  error  term 
in  Eq  (3-6)  indicates  that  the  second  derivative  is  best  approx- 
imated at  the  grid  points.  The  conclusions  reached  for  the 
first  and  second  derivatives  can  be  extended  to  odd  and  even 
derivatives.  That  is,  it  can  be  shown  that  odd  derivatives 
have  the  smaller  error  when  evaluated  at  half  station,  and  even 
derivatives  are  more  accurately  approximated  at  full  station. 
The  trigonometric  half  station  finite  difference  approximation 
for  the  first  derivative  is  derived  in  Appendix  C. 

Appl ication  to  Equilibrium  Equati on 

The  trigonometric  and  conventional  finite  difference 
techniques  will  first  be  applied  to  the  equilibrium  equation 


16 


which  is  repeated  here: 


d^V  . P d^V  _ n 

“• 


(2-1) 


The  derivatives  in  Eq  (2-1)  are  replaced  by  finite  difference 
expressions  which  are  developed  at  each  of  the  internal  ref- 
erence points,  1 through  N,  placed  along  the  beam  as  shown  in 
Fig.  3.3. 


Fig.  3.3  Reference  Points 


This  substitution  in  Eq  (2-1)  results  in  a set  of  N equations, 
one  for  each  reference  point.  These  equations  take  the  form 
of 

7 (Vi-2-''  V,-4 

+ (''^_^-2V.+V.^^)  = 0,  i = 1,2,...,N  (3-36) 

where  Vj^  represents  the  vertical  displacement  at  node  point  k. 
The  boundary  conditions  can  be  applied  by  evaluating  the 
appropriate  finite  difference  expression  at  the  boundary.  For 
example,  a beam  pinned  at  the  left  end  has  boundary  conditions 
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of  V(0)  = 0 and  V"(0)  = 0.  In  terms  of  finite  difference  ex- 
pressions these  two  boundary  conditions  can  be  expressed  as 


or 

and 


V"(0)  = V_^-2Vjj+V^  = 0 

(3-37) 

V-l  = 2V^-Vi 

(3-38) 

V,  = 0 

(3-39) 

The  last  two  equations  relate  Vq  and  v_^  to  the  internal 
displacements  and  provide  the  additional  information  needed  to 
solve  the  system  of  equations  represented  by  Eq  (3-36). 

Since  pinned  and  clamped  beams  are  restrained  in  the  ver- 
tical direction  at  the  boundary,  the  equilibrium  equation  is 
only  applied  at  internal  node  points.  Guided  and  free  beams, 
on  the  other  hand,  may  experience  vertical  deflections  at  the 
boundary.  Therefore,  the  equilibrium  equation  is  applied  at 
the  boundary,  and  is  treated  as  an  unknown.  This  procedure 
introduces  one  additional  equation  and  unknown  in  Eq  (3-36) 
for  beams  which  are  free  or  guided  at  one  end. 

The  N equations  (N  + 1 equations  for  beams  which  are  free 
or  guided  at  the  boundary)  making  up  Eq  (3-36)  can  be  repre- 
sented in  matrix  form 

(3-40) 

where  [C]  is  a function  of  P.  Since  Eq  (3-40)  represents  a 
set  of  homogeneous  algebraic  equations,  it  is  apparent  from 
Cramer's  rule  that  the  determinant  of  the  C matrix  must  be  zero 
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for  a non-trivial  solution.  Therefore,  the  critical,  or  buck- 
ling, force  is  that  value  of  the  axial  force,  P,  which  produces 
a singular  matrix.  A simple  numerical  algorithm  which  can  be 
used  to  find  this  value  of  the  axial  force  involves  an  iterative 
procedure  in  which  the  determinant  is  evaluated  for  successively 
larger  values  of  P.  When  the  determinant  changes  sign,  the 
step  size  is  decreased  and  the  process  repeated  until  the  crit- 
ical force  is  located  to  the  desired  degree  of  accuracy.  This 
eigenvalue  problem  can  also  be  solved  by  more  sophisticated 
techniques  involving  similarity  and  orthogonal  transformations. 
Examples  of  these  methods  include  the  Jacobi  method,  the  Power 
method.  Householder's  method,  and  the  LR  and  QR  algorithms 
[21,25,26].  However,  the  iterative  algorithm  described  above 
has  the  characteristic  of  simplicity  and,  as  such,  lends 
insight  to  the  overall  operation. 

As  mentioned  in  earlier  sections  of  this  report,  an 
initially  flat  beam  with  an  axially  applied  force  will  remain 
flat  for  values  of  the  axial  force  below  the  critical  level. 

When  the  critical  force  is  reached,  the  beam  will  buckle  with 
a characteristic  mode  shape.  This  critical  force  is  the 
lowest  eigenvalue  from  Eq  (3-40)  and  the  mode  shape  is  the 
corresponding  eigenvector.  From  a theoretical  standpoint, 
there  are  an  infinite  number  of  eigenvalues.  That  is,  there 
are  many  values  of  P which  satisfy  Eq  (3-40)  for  non-zero  V(x), 
and  each  eigenvalue  produces  a different  theoretical  mode  shape. 
For  example,  the  eigenvalues  or  buckling  loads  for  a beam 
pinned  at  both  ends,  as  determined  by  analytical  solution  of 
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the  equilibrium  equation,  are 


W f 


n2„2rT 

r„  = n = 1,2,3,...  (3-41) 

n 

The  eigenfunctions  corresponding  to  the  first  and  second  eigen- 
values are  depicted  in  Fig.  3.4. 


Fig.  3.4  First  (Top)  and  Second  (Bottom)  Eigenfunctions 
for  Pi nned-Pi nned  Beam 

Since  buckling  occurs  at  the  smallest  value  of  the  buckling 
load,  the  higher  eigenvalues  have  no  physical  significance 
unless  the  beam  is  constrained  at  the  proper  points.  For  pur- 
poses of  analysis,  however,  these  higher  modes  provide  a means 
for  investigating  more  complex  functions.  In  terms  of  the 
iterative  algorithm  discussed  earlier,  the  procedure  for  find- 
ing these  higher  eigenvalues  is  to  simply  search  for  a value 
of  P larger  than  the  first  eigenvalue  which  again  forces  the  C 
matrix  in  Eq  (3-40)  to  be  singular.  This  is  the  second  eigen- 
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value  and  the  corresponding  values  of  v^  make  up  the  second 
eigenvector.  The  process  can  be  repeated  for  third  and  higher 
ei genval ues . 

The  preceding  analysis,  involving  the  computation  of  the 
first  three  eigenvalues  for  a beam  with  a variety  of  boundary 
conditions  was  performed  using  both  trigonometric  and  poly- 
nomial finite  difference  expressions.  Various  values  of  X were 
used  in  an  effort  to  identify  the  optimum  value  for  each  par- 
ticular boundary  condition  and  eigenvalue.  Finally,  the  trig- 
onometric and  conventional  finite  difference  methods  were 
compared  using  the  following  criteria: 

1.  Required  computer  time 

2.  Degree  of  accuracy  (compared  to  analytical  solution) 
for  a given  number  of  node  points 

3.  Ability  to  predict  optimum  value  of  X 


Appl i cati on  to  Virtual  Work  Equation 

A similar  analysis  \;as  conducted  using  the  virtual  work 
concept.  The  virtual  work  equation  is  repeated  here  for  ease 
of  reference: 


EI 


j: 


dx  dx 


dx 


(2-2) 


0 dx"  dx‘ 

Finite  difference  expressions  are  substituted  for  the  first 
and  second  derivatives  in  Eq  (2-2)  resulting  in 

EI  (j^)  (V_-,-2Vjj+VT)(5V.-,-25VQ+6V^)dx 

= P t-.-r)  {!;  {-\VV^)(-6V^t6V^)dx 


(3-42) 
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Notice  that  half  station  finite  difference  approximations  are 
used  for  the  first  derivative  whereas  the  second  derivative  is 
approximated  at  full  station.  By  assuming  that  E and  I are 
equal  to  unity  and  simplifying  the  expression  for  the  inte- 


grands , 

Eq  (3-42) 

can  be  written  as 

(i)  f dx  = P [*’  g dx 

(3-43) 

where 

f 

= (V_^-2V^+V^)(6V_^-26V^+6V^) 

(3-44) 

g 

= (-Vjj+Vt)(-6V^+6V^) 

(3-45) 

Integration  in  Eq  (3-43)  is  performed  using  the  traper.oid  rule 
to  yield 


1 

7 

L 

N 1 

N 

2 

1 f,) 
i = l ’ J 

= P 

^ 9in/2 

1 = 1 

(3-46) 


where  the  subscripts  on  f and  g denote  the  position  along  the 
beam  (Refer  to  Fig.  3-3)  at  which  Eqs  (3-44)  and  (3-45)  are 
evaluated.  For  example,  assume  the  function  f from  Eq  (3-43) 
can  be  represented  by  Fig.  3.5. 


Fig.  3.5  Trapezoid  Rule 
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The  area  under  the  curve  from  x = 0 to  x = .5  is  approximated 
in  Eq  (3-43)  by  yhf  , and  the  area  under  the  curve  from  x = .5 
to  X = 1.5  is  approximated  by  hf ^ . This  process  is  continued 
to  the  end  of  the  beam.  The  integration  on  the  right  side  of 
Eq  (3-46)  is  performed  in  a similar  manner  with  the  exception 
that  the  function  is  evaluated  at  half  station  points.  The 
summations  in  Eq  (3-46)  can  be  expanded  and  the  coefficients 
of  specific  virtual  di spl acements  , 6V^ , can  be  combined  to  pro- 
duce an  equation  of  the  form 


ta„(P.y 


-1  •''o 


a_i  (P>V_i  »V^,. . . 


or 


N+2 

I = 0 

i = -l  ^ ^ 


(3-47a) 

(3-47b) 


where  are  a function  of  the  actual  displacements,  V(x),  and 
the  axial  force,  P.  It  is  apparent  that  the  internal  virtual 
displacements  in  Eq  (3-47)  are  arbitrary.  However,  the  virtual 
work  concept  as  applied  to  boundary  conditions  states  that  vir- 
tual displacements  at  the  boundary  are  compatible  with  the  real 
boundary  displacements.  For  example,  a beam  pinned  at  the  left 
boundary  must  satisfy  the  conditions  that  6V(0)  = 0 and 
6V"(0)  = 0.  In  terms  of  finite  difference  expressions,  it  is 
necessary  to  insure  that 


0 


-6V 


1 


and 
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(3-48) 

(3-49) 


These  additional  equations  provide  a means  of  expressing  the 
external  virtual  displacements  and  boundary  displacements  in 
terms  of  the  internal  displacements.  Thus,  Eq  (3-47)  can  be 
expressed  as 

N 

I b.6V.  = 0 (3-50) 

i = l ^ ’ 

where  b^  are  also  functions  of  v(x)  and  P.  Since  the  internal 
virtual  displacements  are  arbitrary  in  Eq  (3-50),  the  coeff- 
icients, b^. , must  each  be  zero  for  the  equality  to  hold. 
Therefore,  the  virtual  work  equation  can  be  reduced  to  a set 
of  N equations  which  is  a function  of  the  actual  displacements 
V,  and  the  axial  force,  P.  That  is, 

b^.  = 0 i = 1 ,2,. . . ,N.  (3-51) 

The  boundary  conditions  can  again  be  used  to  relate  the  actual 
external  displacements  in  b^.  to  the  internal  displacements. 

By  factoring  out  the  actual  displacements,  Eq  (3-51)  can  be 
represented  in  matrix  form: 


where  the  elements  of  the  C matrix  are  a function  only  of  P. 

Equation  (3-52)  is  identical  in  form  to  the  matrix  equa- 
tion which  resulted  from  employing  the  equilibrium  approach. 
Thus,  the  same  iterative  technique  can  be  used  here  to  solve 
for  the  critical  force.  It  should  be  emphasized  that  this 
iterative  algorithm  is  referenced  because  of  its  simplicity. 
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The  more  sophisticated  techniques  mentioned  in  the  preceding 
paragraph  are  much  more  efficient,  and  some  of  these  techniques 
have  the  additional  capability  of  computing  eigenvectors  as 
well.  The  methods  are  discussed  at  length  in  numerous  texts 
[21  ,25-27]. 

Mode  Shape  Analysis 

Stein  and  Housner  developed  the  trigonometric  approach  to 
the  finite  difference  calculus  in  the  expectation  that  the 
Fourier  series  would  provide  a closer  approximation  to  the 
actual  displacement  function  than  an  equal  number  of  terms  of 
the  Taylor  series  [1].  To  test  this  hypothesis  and  provide  a 
broader  basis  for  understanding  this  approach,  a modal,  or 
eigenvector,  analysis  was  conducted.  The  purpose  of  this  anal- 
ysis was  to  determine  the  buckled  mode  shape  as  predicted  by 
both  the  trigonometric  and  conventional  approaches.  These  pre- 
dicted mode  shapes  can  be  compared  to  the  theoretical  as 
determined  by  analytical  solution  of  the  buckling  problem. 

The  predicted  mode  shapes  are  represented  by  the  eigen- 
vectors of  Eq  (3-40)  or  (3-52)  which  is  repeated  here: 

= 0 (3-40) 

Recall  that  the  eigenvalue,  P,  is  embedded  in  C.  There  are 
numerous  techniques  available  for  computing  eigenvectors  once 
the  eigenvalue  is  found.  The  most  natural  method  is  to  delete 
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one  row  of  the  C matrix  and  arbitrarily  set  one  element  of  , 
such  as  V-i,  equal  to  unity.  This  reduces  the  order  of  the  C 
matrix  by  one  and  results  in  a set  of  ncn-homogeneous  equations 
for  v.'hich  , i = 2,3,. ..,N,  can  be  computed.  However,  this 
technique  is  known  to  be  inaccurate,  and  the  result  may  be  de- 
pendent on  which  of  the  rows  in  the  matrix  is  deleted  [26]. 

More  promising  algorithms  for  computing  eigenvectors  include 
the  inverse  iteration  method,  the  adjoint  method,  and  the  Jaco- 
bi method  [21,26-28].  These  methods  are  discussed  at  length 
in  the  literature  and  will  not  be  repeated  here. 

To  illustrate  the  nature  of  the  eigenvalue  problem  in 
more  depth,  the  critical  force  can  be  factored  from  the  C 
matrix  in  Eq  (3-40)  to  produce  an  equation  of  the  form 


(3-53) 


This  is  the  representation  for  the  general  eigenvalue  problem. 
Prepared  subroutines  are  available  on  the  International  Mathe- 
matical and  Statistical  Libraries  (IMSL)  for  computing  the 
eigenvectors  of  Eq  (3-53).  These  subroutines  use  similarity 
and  orthogonal  transformations  described  by  Hornbeck  [21]  to 
condition  the  matrices. 

Free  and  guided  beams  introduce  special  problems  in  the 
computation  of  eigenvectors  due  to  the  rigid  body  motion  at 
the  boundary.  This  motion  allows  the  eigenvector  (which  rep- 
resents the  displacement  of  points  along  the  beam)  to  translate 
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arbitrarily  in  the  vertical  direction.  To  specify  the  position 
of  the  reference  line  defining  the  initial  position  of  the 
eigenvector,  the  vertical  displacement  at  the  left  boundary  can 


be  arbitrarily  set  equal  to  zero.  As  a result,  the  computed 
eigenvector  represents  vertical  displacements  from  the  zero 
reference  line,  and  the  rigid  body  motion  is  effectively  elim- 
inated from  the  problem.  If  the  rigid  body  motion  is  allowed 
to  occur,  the  eigenvector  is  produced  in  a form  not  normally 
associated  with  eigenvector  analysis.  That  is,  the  entire 
eigenvector  is  shifted  by  an  amount  proportional  to  the  rigid 
body  motion.  This  translation,  coupled  with  the  normal  scalar 
multiplication  of  eigenvector  elcruents,  produces  a mode  shape 
which  appears  to  conflict  with  the  theoretical  mode  shape. 
Although  it  is  true  that  tfio  computed  and  theoretical  eigenvec- 
tors are  not  the  same,  it  can  be  shown  that  the  correspondi ng 
modal  shapes  are  proportional.  To  avoid  this  complication, 
the  position  of  the  reference  line  can  be  specified  as  outlined 
above . 

Two  independent  techniques  v'ere  used  to  calculate  tiie 
eigenvectors  of  Eq  (3-53).  The  first  of  these  techniques,  the 
inverse  iteration  method,  is  explained  in  detail  by  Franklin 
[26].  Briefly,  this  method  replaces  the  zero  vector  on  tiie 
right  side  of  Eq  (3-40)  by  an  arbitrary  vector,  b.  The 
resulting  set  of  non-homogeneous  equations  can  be  solved  for 
v.j  , i - 1,2,...,N,  the  first  estimate  of  the  eigenvector  which 
replaces  b on  the  right  side  of  Eq  (3-40).  This  process  can 
be  continued  until  convergenc  e IS  a c hieved.  The  second  technique 
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for  computing  the  eigenvectors  involved  the  use  of  a subroutine 
from  IMSL.  The  algorithm  used  by  this  subroutine  was  developed 
by  Moler  and  Stewart  [29].  The  A matrix  in  Eq  (3-53)  is  re- 
duced to  upper  Hessenberg  form  (a  matrix  in  which  the  elements 
a^j  = 0,  for  j < i - 2)  and  the  B matrix  is  reduced  to  upper 
triangular  form.  The  eigenvalues  are  then  computed  from  the 
characteristic  equation  by  a sequence  of  sophisticated  opera- 
tions on  the  elements  of  the  transformed  A and  B matrices. 

The  eigenvectors  are  found  by  using  an  extension  of  the  LR  and 
QR  tri angul ari zations  [21].  These  metliods  are  based  on  an 
iterative  scheme  in  which  the  original  matrices  are  decomposed 
into  the  product  of  lower  triangular  and  upper  triangular 
matrices . 

In  conjunction  with  this  eigenfunction  analysis,  an  inves- 
tigation was  conducted  to  determine  the  accuracy  with  which 
the  truncated  Fourier  and  Taylor  series  are  able  to  predict 
the  actual  boundary  conditions  of  the  buckling  problem.  This 
is  a necessary  prerequisite  for  determining  the  series  approx- 
imation which  most  accurately  represents  the  buckled  function. 
To  accomplish  this  task,  expressions  were  found  for  t!ie  coeff- 
icients of  the  Fourier  series  represented  by  Eq  (3-16): 

it(x-x  ) T^(x-x  ) 

V(x)  = + Tg  sin  + T3  cos  3-I6) 

where  x^  is  the  reference  point.  The  following  expressions 
were  developed  from  a Guass  elimination  scheme: 

f 

I 

I 
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Tl  ■ V,  - ITTolorfT  [V„2-V,-2V,,,cose*2V,cose]  (3-54) 

Tj  " ^^,“5  [-.5V,^2-.5V,+V,^,cose+2V,cos0]  (3-65) 

^3  ' zr^krT  [V(*2-''l-2Vl4l'“*9*2V,cose]  (3-66) 

where  0 = irh/X  and  i represents  the  position  of  the  reference 
point.  Substitution  of  these  expressions  for  the  coefficients 
in  Eq  (3-16)  and  setting  the  reference  point,  x^,  at  node 
point  1 (see  Fig.  3.3)  provides  a means  of  solving  for  V at 
the  boundary  in  terms  of  V at  node  points  1,2,  and  3.  By 
using  the  computed  value  of  the  displacements  at  node  points 
1,2,  and  3,  the  result  obtained  for  gives  an  indication  of 
the  accuracy  with  wtiich  Eq  (3-16)  can  predict  the  particular 
boundary  condition  of  V(0)  = 0.  Boundary  conditions  of 
V'(0)  = 0,  V'‘(0)  = 0,  and  V"  (0)  = 0 can  be  analyzed  in  a 
similar  manner  by  evaluating  the  first,  second,  and  third 
derivatives  of  Eq  (3-16)  at  the  boundary. 

The  same  technique  could  be  used  to  investigate  the 
accuracy  of  the  truncated  Taylor  series.  However,  it  can  be 
shown  that  the  results  obtained  from  the  Fourier  series  equal 
the  results  of  the  Taylor  series  as  X approaches  infinity.  As 
a result,  the  Taylor  series  can  be  analyzed  by  making  X large 
in  Eq  (3-16). 

In  addition  to  looking  at  points  on  the  boundary,  this 
analysis  can  be  extended  to  provide  an  investigation  of  the 
entire  approximating  function.  The  vertical  displacement, 
V(x),  at  arbitrary  points  along  the  beam  can  be  calculated  in 
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terms  of  V^,  Vg,  and  V3  by  using  Eq  (3-16).  In  this  manner  it 
is  possible  to  determine  the  value  of  X which  yields  the  most 
accurate  Fourier  series  expression.  This  concept  will  be  ex- 
panded and  applied  in  the  next  section. 
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IV.  Numerical  Results  and 


Introduction 

The  following  boundary  conditions  were  investigated  using 
both  the  equilibrium  and  virtual  work  approaches: 

1.  Pi nned-Pinned 

2.  Clamped-Pinned 

3.  Cl amped-Cl amped 

4.  Free-Pinned 

5.  Gui ded-P i nned 

6.  Clamped-Free 

7.  Clamped-Guided 

8.  Gui ded-Gui ded 

9.  Free-Free 

10.  Free-Guided 

See  Appendix  D for  a discussion  of  these  boundary  conditions 
and  the  numerical  technique  for  implementation.  Solutions  were 
obtained  for  various  values  of  X ranging  from  .25  to  3.0. 
Without  loss  of  generality,  it  was  assumed  that  the  flexural 
rigidity,  El,  and  the  beam  length,  L,  equal  unity.  The  com- 
puted value  of  the  critical  force  was  compared  to  the  theore- 
tical critical  force,  and  the  percent  error  was  plotted  against 
values  of  X for  each  boundary  condition.  The  mode  shape  anal- 
ysis is  also  presented  in  graphical  form  for  selected  boundary 
conditions.  These  graphs  depict  the  vertical  displacements  at 
buckling  versus  horizontal  points  along  the  beam.  The  boundary 
condition  analysis  was  extended  to  provide  a technique  for 
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calculating  the  Fourier  series  approximating  function.  This 
function  was  plotted  for  various  values  of  A and  arbitrarily 
selected  boundary  conditions. 


Virtual  Work  Approach . 

Calculation  o£  The  critical  force  for  each  boundary 

condition  was  calculated  using  various  values  of  X.  The  per- 
cent error  in  the  computed  value  of  was  plotted  against 
values  of  X as  shown  in  Figures  4.5  through  4.14.  Ten  grid 
points  were  used  for  this  exercise.  Since  the  polynomial,  or 
conventional,  technique  is  independent  of  X,  the  error  result- 
ing from  this  method  appears  as  a horizontal  line.  It  is 
obvious  that  the  various  boundary  conditions  display  similar 
characteristics.  As  a result,  the  discussion  will  be  directed 
toward  tv.-o  of  these  boundary  conditions,  the  pi nned-pi nned 
beam  and  the  free-guided  beam.  The  comments  made  will  be 
applicable  to  the  other  boundary  conditions  as  well. 

The  results  for  the  pinned-pinned  beam  are  shown  in  Fig. 
4.5.  The  dashed  line  in  this  figure  is  used  as  a zero  error 
line  for  ease  of  comparison.  There  is  a range  in  values  of  X 
for  which  the  trigonometric  approach  yields  better  results 
than  the  conventional  approach.  This  range  is  approximately 
.75  < X < ".  It  can  be  seen  that  the  trigonometric  and  conven- 
tional techniques  produce  the  same  P^^  for  large  X as  was 

predicted  in  Section  III.  The  smallest  error  occurs  for  X = L. 

-8 

The  error  in  this  case  is  less  than  1.0  x lo  . (Note  that 
for  this  analysis,  the  beam  length,  L,  was  assumed  equal  to 
1.0).  There  is  a simple  geometric  interpretation  for  this 
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optimum  value  of  The  theoretical  buckled  mode  shape  for 
a pi nned-p i nned  beam  is  shov;n  in  Fig.  4.1. 


Fig.  4.1  Buckled  Mode  Shape  for  Pinned-Pinned  Beam 

The  half  wavelength  of  this  buckled  beam  is  also  equal  to  the 
length  of  the  beam.  Thus,  it  appears  that  the  optimum  value 


of  X corresponds  to  the  half  wavelength  of  the  buckled  mode 


shape.  This  hypothesis  is  supported  by  the  fact  that  X appears 


as  a wavelength  parameter  in  the  original  Fourier  series  of 

: j 

Eq  (3-16)  which  is  repeated  here:  ' 

' « 

tt(x-x  ) tt(x-x  ) I' 

V(x)  = + T2  sin 5^  + T3  cos  (3-16)  ;] 

‘I 

The  results  for  a beam  which  is  free  at  one  end  and  guided  at  i 

the  other  is  shown  in  Fig.  4.6.  The  trigonometric  approach  I 

provides  better  results  than  the  conventional  approach  for 

1.5  < X < «»,  and  the  optimum  value  is  X = 2L.  The  theoretical 

buckled  shape  for  a free-guided  beam  is  depicted  in  Fig.  4.2,  ‘ 

and  the  half  wavelength  is  seen  to  be  equal  to  2L.  Again,  the  ( 

,1 

optimum  value  of  X corresponds  exactly  to  the  half  wavelength  'j 

of  the  buckled  beam.  The  same  conclusion  can  be  drawn  for  all  | 

L 
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F1g.  4.2  Buckled  Mode  Shape  for  Free-Guided  Beam 

boundary  conditions;  and  as  will  be  shov/n  later,  the  same  rule 
applies  to  higher  eigenvalues  as  well.  A summary  of  results 
for  the  virtual  work  approach  with  five  grid  points  is  presen- 
ted in  Table  II  in  Appendix  E.  The  optimum  value  of  X and 
the  range  in  values  of  X for  which  the  trigonometric  approach 
is  superior  to  the  conventional  approach  are  listed  in  the 
first  two  columns.  The  last  two  columns  show  that  X = 1.5 
provides  more  accurate  results  than  the  conventional  approach 
for  all  boundary  conditions. 

Note  in  Fi^gs.  4.5  through  4.14  that  the  calculated  value 
of  the  cr i ti cal ' force  decreases  monoton i ca 1 1 y as  X increases. 
This  can  be  explained  by  examining  the  Fourier  series  approx- 
imating function.  Increasing  the  value  of  X results  in  a 
function  with  a larger  wavelength.  A smaller  force  is  required 
to  produce  a buckled  mode  shape  with  this  increased  wavelength, 
and  the  computed  critical  force  decreases  accordingly. 


Mode  Shape  Analysis.  It  was  anticipated  that  the  value 


of  X yielding  the  best  eigenvalue  would  also  give  the  best 
eigenvector.  Stated  another  way,  it  was  theorized  that  the 
optimum  value  of  X could  be  used  to  predict  the  eigenvalue 
with  a high  degree  of  accuracy  because  the  eigenvector  closely 
represented  the  theoretical  buckled  function.  However,  this 
was  not  entirely  the  case.  Figure  4.15  represents  the  com- 
puted eigenvector  for  a pi nned-pi nned  beam.  A modified  spline 
fitting  technique  [21]  was  used  to  connect  the  elements  of  the 
eigenvector  in  Fig.  4.15.  The  graph  in  the  upper  right  hand 
corner  of  the  figure  represents  the  theoretical  buckled  mode 
shape.  It  can  be  seen  that  the  eigenvectors  calculated  by 

the  polynomial  and  trigonometric  approaches  match  the  theore- 
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tical  mode  shape  within  a tolerance  of  1.0  x 10  . Figure 

4.16  shows  a free-guided  beam  with  identical  results.  The 
same  conclusions  were  reached  for  the  other  boundary  conditions. 
Thus,  the  eigenvector  can  be  calculated  exactly  regardless  of 
the  value  of  X.  This  surprising  result  can  be  partially  ex- 
plained by  looking  at  the  eigenvalue  problem  as  represented  by 
Eq  (3-52): 


0. 


(3-52) 


Recall  that  the  C matrix  is  a function  of  both  X and  the 
eigenvalue,  It  appears  that  the  selection  of  a poor 

value  for  X is  neutralized  by  the  computed  eigenvalue.  That 
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is,  two  distinct  values  of  X result  in  the  calculation  of  two 
different  eigenvalues.  However,  the  C matrix  evaluated  at 
one  value  of  X and  the  corresponding  eigenvalue  differ  only 
by  a scalar  constant  from  the  C matrix  evaluated  at  the  other 
value  of  X and  corresponding  eigenvalue.  This  hypothesis  was 
verified  for  several  specific  cases.  For  example,  the  first 
eigenvalue  of  Eq  (3-52)  with  X = 1.0  and  N = 5 for  a pinned- 
pinned  beam  is  9.869604*,  the  first  eigenvalue  with  X = 2.0  and 
N = 5 is  9.701455.  Evaluation  of  the  C matrix  for  these  two 
cases  results  in  matrices  whose  elements  differ  by  the  con- 
stant factor  of  1.02011.  From  a physical  point  of  view,  the 
equilibrium  and  virtual  work  differential  equations  as  appro- 
ximated by  finite  difference  expressions  (from  which  the 
eigenvalue  problem  was  developed)  require  that  the  buckled 
mode  shape  be  preserved.  It  appears  that  the  shape  oF  the 
approximating  Fourier  series  function  is  close  enough  to  the 
theoretical  shape  to  relate  local  nodal  point  displacements 
with  high  accuracy.  To  accomplish  this  preservation  of  the 
modal  function,  the  eigenvalue  (or  critical  force)  is  adjusted 
such  that  the  eigenvector  remains  constant  regardless  of  the 
value  of  X. 

Fourier  Series  Approximati no  Function . It  is  apparent 
that  the  mode  shape,  or  eigenvector,  analysis  is  unable  to  pro- 
vide an  explanation  for  the  behavior  of  the  finite  difference 
technique  as  formulated  in  this  thesis.  The  next  stop  is  to 
determine  how  closely  the  basic  Fourier  and  Taylor  series 
match  the  actual  boundary  conditions.  The  procedure  used  in 
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this  analysis  was  outlined  in  Section  III.  The  results  (nor- 
malized to  the  beam  length)  for  a beam  pinned  at  both  ends 
are  shown  in  Table  I.  The  boundary  conditions  for  a pinned- 
pinned  beam  are  V(0)  = 0 and  V‘'(0)  = 0.  It  can  be  seen  in  ; 

Table  I that  the  Fourier  series  v.'ith  X = 1.0  matches  these 


Table  I 

Boundary  Condition  Analysis  for  Pinned-Pinned  Beam 


Series 

Representa  ti on 

V(0) 

V"(0) 

Taylor 

.021 

5.30 

Fourier  X = .50 

.Ofil 

13.48 

X = .75 

.016 

3.84 

X =1.0 

-1  5 

5.3x10 

2.4x10"^^ 

X =1.5 

-.012 

-2.89 

X =2.0 

-.016 

-3.94 

boundary  conditions  very  closely.  Recall  that  this  is  the 
optimum  value  of  X for  a pinned-pi nned  beam.  It  can  also  be 
noted  from  Table  I that  the absol ute  val ue  of  the  Taylor  series 
error  is  less  than  the  Fourier  series  error  with  X = .5  but 
more  than  the  Fourier  series  error  with  X = .75.  This  corres- 
ponds exactly  to  the  error  in  calculating  P^^.  The  absolute 
value  of  the  error  using  the  polynomial  approach  is  .68%  com- 
pared to  2.1%  for  the  trigonometric  approach  with  X = .5  and 
.53%  with  X = .75.  Thus,  there  is  a direct  correlation 
between  the  accuracy  in  calculating  P^^  and  the  accuracy  with 
which  the  assumed  function  satisfies  the  boundary  conditions. 

This  boundary  condition  analysis  can  be  extended  to 
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provide  an  investigation  of  the  entire  function  represented 
by  the  Fourier  series.  From  Eqs  (3-54)  through  (3-56),  it  was 
shown  that  the  coefficients  of  the  Fourier  series  can  be  ex- 
pressed in  terms  of  X and  the  displacements  at  three  neighbor- 
ing points  along  the  beam.  If  these  coefficients  are  evaluated 
for  a particular  value  of  X and  substituted  in  the  Fourier 
series,  a functional  relationship  is  obtained  in  which  x and 
v(x)  are  the  only  variables: 

tt(x-x^)  it(x-x  ) 

V(x)  = sin + T3  cos  (3-16) 

where  x^  represents  the  point  about  which  the  Fourier  expan- 
sion is  made.  Note  that  V(x)  is  dependent  on  the  values  of 
the  three  computed  eigenvector  displacements  as  outlined  in 
Section  III  which  ore  used  to  compute  the  coefficients.  Using 
this  technique,  it  is  possible  to  determine  the  function  being 
represented  by  the  Fourier  series.  It  is  important  to  dis- 
tinguish between  V(x)  as  represented  by  Eq  (3-16)  and  the  V(x) 
which  was  found  from  interpolation  of  the  eigenvector  compon- 
ents. The  former  is  the  original  assumed  function  from  which 
finite  difference  expressions  were  developed.  The  latter 
represents  the  results  of  an  eigenvector  analysis.  The  primary 
difference  is  that  the  first  is  an  assumed,  or  approximating 
function,  whereas  the  second  is  a calculated  function  based 
upon  this  assumed  expression. 

It  has  already  been  shown  that  the  eigenvector  of  the 
finite  difference  matrix  equation  is  aligned  perfectly  with 
the  theoretical  eigenvector  regardless  of  the  value  of  X. 
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The  next  step  is  to  determine  if  a relationship  exists  between 
A and  the  approximating  function  of  Eg  (3-16).  It  would  be 
reasonable  to  anticipate  that  the  assumed  function  ir.corpor- 
I ating  the  value  of  X which  yields  the  best  approximation  to 

' the  actual  buckled  shape  will  correspond  to  tfie  optimum  value 

of  X for  computing  the  critical  force.  To  test  this  hypothesis, 

Eq  (3-16)  was  evaluated  for  a beam  pinned  at  both  ends.  The 
resulting  functions  for  two  values  of  X are  shown  with  the 
corresponding  Taylor  series  expansion  in  Fig.  4.17.  The 
theoretical  buckled  function  i s exactly  reproduced  and  boundary 
conditions  are  pei^fectly  matched  when  a value  of  1.0  is  used 
for  X.  The  results  for  a beam  which  is  free  at  one  end  and 
guided  at  the  other  are  sliov/n  in  Fig.  4.18.  In  this  case, 
the  function  corresponding  to  X = 2.0  (the  optimum  value  cf  X 
for  a free-guided  beam)  matches  the  theoretical  buckled  mode 
shape  and  boundary  conditions.  Identical  results  were  obtained 
for  the  remaining  boundary  conditions  v/hich  are  shov.n  in  Figs. 

4.19  through  4.26. 

It  follows  from  this  analysis  that  the  accuracy  of  the 
trigonometric  and  conventional  finite  difference  methods  is 
directly  related  to  the  accuracy  with  which  the  approximating 
function  represents  the  buckled  mode  shape  and  boundary  condi- 
tions. 

Higher  Order  Eigenvectors.  As  mentioned  previously, 
second  and  higiter  order  eigenvalues  of  the  buckling  equation 
have  no  physical  significance.  From  a theoretical  standpoint, 
iiowever,  these  eigenvalues  are  valuable  since  the  corresponding 
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characteristic  function  tends  to  be  more  complex.  The  results 
of  computing  the  second  eigenvalue  for  a beam  pinned  at  each 
end  are  shown  in  Fig.  4.27.  Comparison  of  this  graph  with 
Fig.  4.5,  which  depicts  the  corresponding  first  eigenvalue, 
indicates  that  the  error  tends  to  be  higher  for  the  second 
eigenvalue.  Note  that  the  optimum  value  of  A is  now  .5.  The 
second  buckled  mode  for  a beam  pinned  at  both  ends  is  shown 
in  Fig.  4.3  and  shows  that  A = .5  still  corresponds  to  the 
half  wavelength.  The  results  of  the  eigenvector  computation 
is  shown  in  Fig.  4.28. 

The  results  of  calculating  the  third  eigenvalue  for  a 
pi nned-pi nned  beam  are  shown  in  Fig.  4.29.  The  optimum  value 


Fig.  4.3  Second  Buckled  Mode  for  a Pi nned-Pinned  Beam 

of  A is  .33  which  again  corresponds  to  the  half  wavelength  as 
indicated  by  Fig.  4.4.  The  eigenvector  analysis  is  displayed 
in  Fig.  4.30.  The  results  of  calculating  higher  eigenvalues 
for  the  remaining  boundary  conditions  display  identical 
characteristics. 


i 


Fig.  4.4  Third  Buckled  Mode  for  a Pinned-Pinned  Beam 

Degrees  of  Freedom.  The  truncati on  error  of  the  finite 

difference  method  inci'eases  as  the  number  of  grid  points  is 

decreased.  To  quantify  this  change  in  the  truncation  error, 

the  critical  force  for  a pinned-pinned  beam  was  computed  using 

r five  grid  points;  and  the  results  are  shown  in  Fig.  4.31.  As 

expected,  the  error  is  approximately  3.35  times  the  magnitude 

of  the  error  resulting  from  the  use  of  ten  grid  points,  but 

the  optimum  value  of  X remains  at  1.0.  It  is  possible  '.o  carry 

this  analysis  one  step  further.  Recall  that  the  central  finite 

difference  expressions  used  in  this  thesis  have  an  error  term 

2 

which  IS  of  order  h . For  ten  grid  points,  h = .091  and 
h^  = .0083;  for  five  grid  points,  h - .167  and  h^  = .0278.  It 
follows  that  the  error  term  for  five  grid  points  should  be 
approximately  3.36  times  larger  than  the  error  term  for  ten 
grid  points.  This  ratio  was  substantiated  computationally  for 
all  boundary  conditions.  As  a further  example  of  this  charac- 
teristic, Fig.  4.32  depicts  the  results  of  using  five  grid 
points  for  a free-guided  beam.  The  average  ratio  of  the  error 
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terms  for  five  grid  points  compared  to  ten  grid  points  is  3.38. 

It  was  shown  in  the  preceding  section  that  the  truncation  error 

for  trigonometric  finite  difference  approximations  is  the  same 

as  the  error  for  conventional  approximations  when  X approaches 

2 

infinity.  The  error  in  both  cases  is  of  order  h . Although  a 
formal  proof  is  not  presented,  it  appears  from  the  analysis  con- 
ducted above  that  the  error  term  for  trigonometric  approximations 
2 

is  of  order  h for  all  values  of  X.  It  is  also  important  to 
note  that  a very  small  number  of  grid  points  can  be  used  to 
calculate  a highly  accurate  critical  force  if  a judicious  choice 
is  made  for  X.  Of  course,  it  is  also  true  that  a very  large 
error  can  result  from  a poor  choice  for  x.  These  conclusions 
were  found  to  be  equally  valid  for  all  boundary  conditions 
investigated  in  this  thesis.  A summary  of  results  with  N = 5 
is  presented  in  Table  II  of  Appendix  E. 

Full  Station  vs  Half  Station.  It  was  sho\/n  in  Section 
III  that  odd  derivatives  are  more  accurately  represented  by 
half  station  finite  difference  approximations.  As  a result, 
the  virtual  work  approach  was  implemented  by  replacing  the 
first  derivative  with  the  finite  difference  approximation  eval- 
uated at  half  station  points.  To  investigate  the  magnitude  of 
accuracy  degradation,  the  virtual  work  approach  was  reformulated 
using  the  full  station  finite  difference  approximation  for  the 
first  derivative.  The  results  of  tliis  trial  for  a pinned- 
pinned  beam  using  ten  grid  points  are  shown  in  Fig.  fl.33.  Com- 
parison with  Fig.  4.5  indicates  that  the  use  of  full  station 
points  to  calculate  the  first  derivative  results  in  an  increase 
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in  the  percent  error  of  approximately  .B%.  The  same  analysis 
for  a free-guided  beam  resulted  in  an  increase  in  the  percent 
error  of  approximately  4%.  In  addition  to  the  decrease  in 
accuracy,  it  should  also  be  noted  that  the  optimum  value  of  X 
has  shifted  from  1.0  to  approximately  .5  for  a pi nned-pi nned 
beam  and  from  2.0  to  a value  less  than  .25  for  a free-guided 
beam.  It  appears  that  the  optimum  value  of  X corresponds  to 
the  half  v.-avelength  of  the  modal  shape  only  v/hen  the  most  accur- 
ate (smallest  truncation  error  term)  finite  difference  repre- 
sentations are  used. 

Computer  Time.  The  amount  of  required  computer  time  vvas 
found  to  be  approximately  the  same  for  both  the  polynomial  and 
trigonometric  approaches.  For  example,  the  polynomial  techni- 
que required  .157  seconds  of  computer  time  to  compute  both  the 
critical  force  and  the  eigenvector  for  a pi nned-pinned  beam 
v/ith  ten  i.ode  points;  the  trigonometric  technique  with  X = l.O 
required  .147  seconds.  The  difference  between  these  two  mea- 
surements is  well  within  the  accuracy  of  the  computer  clock 
and,  as  such,  can  be  considered  negligible.  The  same  conclusion 
was  reached  for  all  boundary  conditions.  The  results  for  five 
node  points  were  similar.  The  polynomial  technique  required 
.025  seconds;  the  trigonometric  technique  with  X - 1.0  required 
.024  seconds. 


Vi rtiia  1 Work  and  G a 1 e r k i n Approaches  . It  is  interesting 
to  note  that  the  virtual  work  approach,  as  incorporated  in  this 
thesis  is  similar  to  the  Galerkin  approach.  Both  methods  are 
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based  on  the  principle  of  minimum  potential  energy.  The  form 
of  the  Galerkin  method  as  incorporated  herein  expresses  the 
potential  energy  of  the  system  by  applying  a virtual  displace- 
ment to  the  differential  equation  which  describes  the  equili- 
brium of  all  forces  acting  on  the  system.  This  expression  is 

•L 

[L(V)-P]6V  dx  = 0 (4-1) 

where  L denotes  a differential  operator  acting  on  the  displace- 
ment [30].  The  unknown  displacements  are  approximated  by  a 
series  function  which  satisfies  the  prescribed  boundary  con- 
ditions: 

M 

V(x)  = I a.(!).(x)  (4-2) 

i = l • 

where  are  the  undetermined  coefficients  and  1>^(x)  repre- 
sent continuous  functions.  The  expression  for  the  virtual 
d i s pi  a cement  i s then 

M 

6V(x)  = I 6a. 4). (x).  (4-3) 

1 = 1 ^ ^ 

Substitution  of  Eqs  (4-3)  and  (4-2)  in  Eq  (4-1),  expansion  of 
the  variational  series  expression,  and  noting  that  the  varia- 
tions of  the  coefficients  are  arbitrary  result  in  the  following 
system  of  equations: 

4.  . dx  = 0 j = 1 ,2,3,.  . . ,M.  (4-4) 

J 

Integration  of  Eq  (4-4)  yields  a final  system  of  equations 

which  can  be  expressed  as 
M 

I f.  :(P)a.  = 0 i = l,2....,?i  (4-5) 

j = l J 


L 

^0 


M 

L(  I 
1 = 1 
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in  which  f..  denotes  a functional  relationship  in  terms  of 

* J 

the  external  force.  Since  the  coefficients,  a^,  are  arbitrary, 
Eq  (4-5)  represents  an  eigenvalue  problem  from  which  the  crit- 
ical force  can  be  computed. 

Although  the  Galerkin  and  virtual  work  approaches  are  not 
identical,  the  similarities  are  obvious.  In  each  case,  the 
energy  equation  can  be  derived  from  equilibrium  expressions, 
and  the  displacement  function  is  represented  by  a series  app- 
roximation. The  Galerkin  approach  results  in  M homogeneous 
equations  with  M arbitrary  coefficients  where  M is  the  number 
of  terms  in  the  series  approximation,  and  herein  lies  the  pri- 
mary difference  between  the  two  methods.  The  virtual  work 
approacli  outlined  in  this  thesis  develops  N equations  with  K 
arbitrary  coefficients  where  M is  the  number  of  grid  poiits. 
These  N equations  are  represented  by  Eq  (3-52)  which  can  be 
expressed  as 

N 

I gii(P)Vi  =0  i = 1,2,...,N.  (4-6) 

j=l  ^ 

In  this  case  the  arbitrary  variables  (or  eigenvector  compon- 
ents) correspond  to  the  displacement  function.  For  the 
Galerkin  method,  the  coefficients  of  the  approximating  func- 
tion are  the  components  of  the  eigenvector  as  seen  in  Eq  (4-5). 
Since  these  coefficients  determine  the  displacement  function, 
v(x),  it  is  apparent  that  the  Galerkin  and  finite  difference 
approximation  of  virtual  work  methods  are  closely  related. 

The  similarity  becomes  more  pronounced  when  the  general 
concept  behind  the  Galerkin  approach  is  analyzed.  This  con- 
cept is  based  on  the  fact  that  the  error  in  (a^<t)^)  is  minimized 
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for  any  value  of  M if  the  a^.  are  chosen  such  that  Eqs  (4-4) 
are  simultaneously  satisfied  [3].  In  the  virtual  work 
approach,  the  coefficients  of  the  Fourier  series  can  be  ex- 
pressed as  a function  of  the  displacement  values,  v^,  as  shown 
by  Eqs  (3-54)  through  (3-56).  The  v^  are  implicitely  selected 
(through  the  computational  properties  of  the  virtual  work  al- 
gorithm) such  that  the  error  due  to  the  approximation  of  v(x) 
is  minimized.  This  supposition  was  substantiated  numerically 
by  the  analysis  of  the  Fourier  series  approximating  function 
presented  in  preceding  paragraphs.  It  was  also  shown  that 
full  station  approximations  for  the  first  derivative  result 
in  larger  errors  in  the  approximation  of  V(x)  and,  consequently, 
larger  errors  in  the  calculation  of  P^^,. 

Equi 1 i bri urn  Approach 

Calculation  of  P ...  Many  of  the  comments  made  about  the 

— 

virtual  work  approach  will  also  be  applicable  to  the  equil- 
ibrium approach.  The  results  of  computing  P^^  for  the  various 
boundary  conditions  are  shown  in  Figs.  4.34  through  4.43. 

These  figures  have  similar  characteristics,  and  the  pinned- 
pinned  and  free-guided  beams  were  again  singled  out  for  dis- 
cussion . 

The  results  using  the  equilibrium  approach  for  a pinned- 
pinned  beam  are  shown  in  Fig.  4.34.  The  trigonometric 
approach  is  more  accurate  than  the  conventional  approach  for 
values  of  X in  the  range  .9  < X < 1.1  and  2.0  < X < <».  As  was 
true  for  the  virtual  work  approach,  the  optimum  value  of  X 
corresponds  to  the  half  wavelength  of  the  buckled  mode  shape. 


For  the  pinned-pinned  beam,  this  occurs  at  X = l.C;  and  the 
error  is  less  than  1.0  x 10  . There  is  a second  value  of  ^ 
for  which  the  error  in  calculating  approaches  zero.  This 

value  for  a pinned-pinned  beam  is  approximately  2.8.  The 
mathematical  explanation  for  this  second  optimal  value  will  be 
given  later. 

The  results  for  a free-guided  beam  are  shown  in  Fig.  4.35. 
There  is  again  a range  in  values  of  X for  which  the  trigono- 
metric approach  gives  more  accurate  results  than  the  conven- 
tional approach,  and  the  optimum  values  of  X are  2.0  and  5.25. 

The  first  of  these  values  corresponds  to  the  half  wavelength 
for  this  boundary  condition  and  provides  accuracy  within 

1.0  >c  10“^. 

A summary  of  results  for  the  equilibrium  approach  with 
five  grid  points  is  presented  in  Table  III  of  Appendix  F.  The 
first  two  columns  list  the  optimum  values  of  >.  and  the  range 
in  values  of  X for  which  the  trigonometric  approach  is  superior 
to  the  conventional  approach.  The  last  two  columns  show  that 
X = 3.75  provides  more  accurate  results  than  the  conventional 
approach  for  all  boundary  conditions. 

It  is  important  to  note  that  the  equilibrium  and  virtual 
work  methods  produce  the  same  value  for  when  conventional 
finite  difference  approximations  are  used.  For  example,  the 
error  in  computing  for  a pinned-pinned  beam  using  the 
equilibrium  method  and  conventional  finite  difference  approxi- 
mations differs  from  the  error  in  computing  P^^  using  the  virtual 
Work  method  by  less  tfian  1.0  x lo"^^ 


This  observation  attests 


to  the  basic  similarity  between  the  two  methods.  The  methods  do 


not  produce  the  same  however,  when  trigonometric  finite 

difference  approximations  are  used.  In  fact,  the  equilibrium 
approach  is  much  more  sensitive  to  errors  in  the  estimate  of  X 
than  the  virtual  work  approach.  For  example,  the  virtual  work 
approach  produces  an  error  of  1.5%  in  the  calculation  of  for 
X = .75  The  equilibrium  approach  produces  an  error  of  7.^% 
for  the  same  X value. 


Mode  Shape  Analysis.  As  was  found  to  be  true  for  the  vir- 
tual work  approach,  the  mode  shape  analysis  was  not  useful  in 
explaining  the  degree  of  accuracy  obtained  from  the  equilib'^ium 
approach  for  different  values  of  X.  The  results  of  this  anal- 
ysis for  a pinned-pinned  beam  are  shewn  in  Fig.  4.44.  It  can 
be  seen  that  the  computed  eigenvector  matches  the  theoretical 
buckled  function  perfectly  for  the  values  of  X presented.  The 
same  results  are  obtained  for  a free-guided  beam  as  shown  in 
Fig.  4.45.  A possible  explanation  for  this  behavior  is  the 
seme  as  that  proposed  for  the  energy  approach. 


Fourier  Series  Approximating  Function.  In  order  to  inves- 
tigate the  Fourier  series  used  by  the  equilibrium  approach,  it 


is  necessary  to  solve  for  five  coefficients: 


V(x)  = T,  + T,  sin 


tt(x-x„) 


+ Tj  cos 


■'t(x-x„) 


+ T^  sin  2 


■rT(  x-x,  ) 


+ Tg  cos  2 


^(x-x  ) 

* r\  • 


(4-7) 


This  is  necessary  since  the  finite  difference  approximation 


for  the  fourth  derivative  in  the  equilibrium  equation  was 
derived  from  Eq  (4-7).  Recall  that  only  three  terms  of  the 
Fourier  series  were  used  in  the  development  of  the  virtual 
work  approach.  Expressions  for  the  coefficients  in  Eq  (4-7) 
were  developed  in  the  same  manner  as  the  development  of  co- 
efficients for  the  three  term  Fourier  series  as  used  in  the 
virtual  work  method.  The  functions  obtained  from  evaluating 
Eq  (4-7)  for  the  various  boundary  conditions  are  shown  in 
Figs.  4.46  through  4.55.  The  results  for  a pi nned-pi nned  beam, 
which  is  depicted  in  Fig.  4.46,  is  characteristic  of  the 
remaining  boundary  conditions  and  will  be  used  as  an  example. 

As  expected,  the  theoretical  buckled  mode  shape  is  matched 
perfectly  by  the  approximating  function  when  evaluated  at 
X = 1.0.  This  value,  of  course,  corresponds  to  the  half  wave- 
length. For  this  particular  function,  T^  and  Tg  were  numer- 
ically determined  to  be  zero.  Although  not  shown  pictorially, 
it  was  found  that  Eq  (4-7)  matches  the  buckled  mode  shape  when 
X = 2 as  well  as  when  X = 1.  The  values  of  ^3 

observed  to  be  zero  for  X = 2.  With  this  information  available, 
it  i s apparent  from  Eq  (4-7)  that  the  expression  for  V (x)  is 
the  same  for  X = 1 and  X = 2.  At  first  glance,  it  would 
appear  that  these  values  of  X should  give  the  most  accurate 
Per"  However,  reference  to  Fig.  4.34  shows  that  X = 1 and 
X = 2.8  are  the  optimizing  values  of  X.  This  apparent  contra- 
diction can  be  partially  explained  by  recalling  that  the 
equilibrium  equation  contains  both  the  fourth  derivative  and 
the  second  derivative.  The  five  term  representation  of  the 


49 


r 


Fourier  series  (used  to  derive  the  fourth  derivative)  matches 
the  buckling  mode  shape  for  X values  of  one  or  two,  whereas 
the  three  term  representation  (used  to  derive  the  second 
derivative)  matches  only  at  A equal  to  one.  It  is  obvious 
that  X = 1 should  be  an  optimum  value  for  calculating  by 
the  equilibrium  approach.  The  second  optimizing  value  of  X is 
the  one  that  provides  the  correct  combination  of  terms  from 
both  the  three  term  and  five  term  Fourier  series.  For  a beam 
pinned  at  both  ends,  this  combination  appears  to  occur  at 
X = 2.8. 

The  important  point  is  that  there  are  two  va'-ues  of  X 
from  which  the  exact  value  of  can  be  obtained  using  the 
equilibrium  approach.  The  first  value  of  X corresponds  to  the  j 

half  wavelength  of  the  buckled  mode  shape.  The  second  optimi- 
zing X is  more  difficult  to  predict  since  the  value  which 
provides  a very  close  approximation  for  the  fourth  derivative 
appears  to  result  in  an  extremely  poor  approximation  for  the 
second  derivative.  The  end  result  is  an  optimum  value  which 
is  larger  than  both  of  the  individual  candidates. 

The  same  conclusions  were  reached  for  each  of  the  remain- 
ing boundary  conditions.  For  example.  Fig.  4.47  illustrates  ■ 

the  approximating  Fourier  series  for  a beam  free  at  one  end 
and  guided  at  the  other.  The  theoretical  mode  shape  is  repro- 
duced exactly  for  X = 2.0  which  correspond’'  to  the  half  wave- 
length for  this  boundary  condition. 

Higher  Order  Eigenvectors.  The  results  of  computing  the 
second  eigenvalue  for  a beam  pinned  at  both  ends  are  shown  in 


50 


Fig.  4.56.  It  is  apparent  that  the  error  in  calculating  the 
second  eigenvalue  is  higher  than  the  error  produced  from  cal- 
culating the  first  eigenvalue.  The  optimum  value  of  X is  .5 
which  corresponds  to  the  half  wavelength  of  the  buckled  mode. 
There  is  a second  optimum  value  at  X = 1.0.  The  corresponding 
eigenvector  is  represented  by  Fig.  4.57.  This  relatively  com- 
plex eigenfunction  is  reproduced  exactly  (standard  deviation 
of  errors  is  less  than  1.0  x 10"^^)  regardless  of  the  value  of 
X . The  conclusions  reached  for  the  second  eigenvalue  extend 
naturally  to  the  third  eigenvalue  which  is  depicted  in  Figs. 

4.58  and  4.59.  The  error  is  higher  for  each  value  of  X than 
for  the  first  and  second  eigenvalues,  and  the  optimum  value  of 
X is  .33  which  again  corresponds  to  the  buckled  half  wavelength. 

Degrees  of  Freedom.  The  results  of  using  five,  rather 
than  ten,  grid  points  for  a beam  pinned  at  each  end  are  shov/n 
in  Fig.  4.60.  As  expected,  the  error  is  approximately  3.34 
times  higher  for  each  value  of  X;  but  the  optimum  value  of  x 
is  still  1.0.  (For  additional  boundary  conditions  incorpor- 
ating N'  = 5,  refer  to  Table  III  of  Appendix  E).  The  computed 
eigenvector  matches  the  theoretical  buckled  function  with  the 

_Q 

standard  deviation  of  errors  less  than  1.0  x lo  . Similar 
results  were  obtained  for  the  other  boundary  conditions. 

Computer  Time.  The  difference  between  required  computer 
time  for  the  trigonometric  approach  and  computer  time  for  the 
conventional  approach  was  found  to  be  insignificant.  For 
example,  tfie  conventional  technique  required  .152  seconds  to 
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compute  the  critical  force  and  eigenvector  for  a pi nned-pi nned 
beam  with  ten  grid  points.  The  trigonometric  technique  with 
X = 1.0  required  .149  seconds.  The  conventional  technique 
required  .026  seconds  with  five  grid  points  compared  to  .025 
for  the  trigonometric  approach  with  X = 1.0. 
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Fig.  4.27  Second  Eigenvalue  for  Pi nned-Pinned  Beam  - N 

= 10  (Virtual  Work) 
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= 10  (Equilibrium) 
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V.  Conclusions 


1 
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This  thesis  compares  the  results  of  using  the  trigonomet- 
ric and  conventional  approaches  to  the  finite  difference 
calculus  to  solve  the  equilibrium  and  virtual  work  equations. 

A wide  range  of  boundary  conditions  were  investigated.  Various 
values  of  X were  used  in  the  trigonometric  approach  to  deter- 
mine the  optimum  value  as  well  as  to  determine  the  range  over 
which  the  trigonometric  approach  gives  more  accurate  approxi- 
mations than  the  conventional  approach.  In  addition,  an  in- 
depth  search  was  conducted  to  provide  plausible  explanations 
for  the  superiority  of  one  method  over  the  other  and  one  value 
of  X over  other  values.  Finally,  the  effect  of  decreasing  the 
number  of  grid  points  and  the  use  of  full  station  approximations 
in  the  virtual  work  equation  were  investigated. 

The  virtual  work  method  was  found  to  be  an  efficient  and 
simple  approach  and  provided  excellent  results  for  both  the 
trigonometric  and  conventional  techniques  with  as  few  as  five 
grid  points.  As  predicted  from  theory,  computational  data  re- 
vealed that  the  magnitude  of  error  in  computing  varies 
directly  with  the  square  of  the  grid  size.  The  variable  input 
parameter,  X,  has  the  effect  of  adjusting  the  wavelength  of 
the  Fourier  series  approximating  function,  and  the  optimum  value 
of  X corresponds  to  the  half  wavelength  of  the  buckled  mode 
shape  for  each  boundary  condition.  There  is  a range  in  the 
values  of  X for  which  the  trigonometric  approach  is  superior. 
This  range  extends  from  approximately  25%  below  the  optimum 
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value  to  infinity.  Thus,  a large  value  of  A is  guaranteed  to 
provide  more  accurate  results  than  the  conventional  approach. 

Of  course,  if  X is  chosen  to  be  too  large,  the  error  from  the 
conventional  and  trigonometric  techniques  approach  the  same 
value;  and  the  benefit  of  using  the  trigonometric  technique  is 
lost.  It  was  shown  in  Appendix  E that  X = 1.5L  produces  satis- 
factory results  for  all  boundary  conditions.  A potential  ex- 
planation for  the  superiority  of  one  method  over  the  other  was 
found  from  an  investigation  of  the  Fourier  and  Taylor  series 
approximating  functions.  The  X value  which  yields  the  closest 
series  approximation  to  the  theoretical  displacement  function 
corresponds  to  the  optimum  value  of  X.  That  is,  the  key  to 
calculating  an  accurate  estimate  of  the  critical  force  is  to 
supply  an  approximating  function  which  very  closely  reproduces 
the  buckled  mode  shape. 

Finally,  the  similarities  between  the  virtual  work  tech- 
nique as  employed  in  this  thesis  and  the  Galerkin  approach 
were  explored.  In  both  cases,  equilibrium  expressions  are  used 
to  derive  potential  energy  relationships;  and  the  displacement 
functions  are  approximated  by  series  expansions.  In  addition, 
there  is  a strong  relationship  between  the  resulting  sets  of 
equations  developed  by  the  two  methods.  Both  methods  attempt 
to  minimize  the  error  in  approximating  the  displacement  func- 
tion. When  the  approximating  function  is  altered  such  that 
this  minimized  error  is  larger,  the  error  in  the  computed  crit- 
ical force  will  increase  proportionally.  This  concept  was 
demonstrated  by  the  use  of  the  full  station  finite  difference 


r 


no 


approximation  for  the  first  derivative.  The  decreased  accuracy 
in  the  approximation  of  v(x)  caused  a significant  increase  in 
the  error  of  the  calculated  value  of  the  critical  force. 

The  equilibrium  approach  was  also  found  to  be  efficient, 
and  excellent  results  were  obtained  using  both  the  trigonometric 
and  conventional  techniques.  However,  it  is  more  difficult  to 
select  an  effective  value  of  X using  this  approach  than  was 
found  to  be  true  for  the  virtual  work  approach.  There  are  two 
optimum  values  of  X,  and  the  first  corresponds  to  the  half 
wavelength  of  the  buckled  mode  shape.  A precise  estimate  of 
the  buckled  wavelength  is  required  in  this  case,  however,  since 
there  is  little  margin  for  error.  The  range  around  this  opti- 
mum value  for  which  the  trigonometric  approach  is  superior  to 
the  conventional  approach  is  very  small,  and  the  error  builds 
rapidly  as  estimates  of  the  optimum  value  worsen.  There  is  a 
large  range  around  the  second  optimum  value  for  which  the  trig- 
onometric approach  is  superior.  This  range  extends  from 
approximately  27%  below  the  optimum  value  to  infinity.  This 
provides  a comfortable  margin  of  error  for  selecting  X.  The 
problem  is  that  there  is  no  known  physical  parameter  from  which 
this  second  optimal  value  can  be  estimated.  It  appears  from 
the  available  data  that  a value  which  is  2.75  times  the  half 
wavelength  provides  a reasonably  close  estimate  in  most  cases, 
but  specific  boundary  conditions  vary  considerably  from  this 
figure.  Despite  the  uncertainty,  it  is  much  safer  to  attempt 
an  estimate  of  the  second  optimal  value  of  X due  to  the  larger 
error  margin.  An  attempt  to  use  the  first  optimal  value  is 


probably  unwise  unless  the  buckled  mode  shape  Is  known  a priori 
with  reasonable  accuracy.  It  was  shown  In  Table  III  of  Appendix 
E that  A = 3.75L  provides  more  accurate  results  than  the  con- 
ventional approach  for  all  boundary  conditions. 

In  comparing  the  results  of  the  virtual  work  and  equil- 
ibrium approaches,  many  similarities  were  noticed  despite  the 
major  conceptual  differences  in  the  derivation  of  these  methods. 
The  interpretation  of  the  wavelength  parameter.  A,  Is  the  same 
In  both  cases  as  already  discussed.  In  addition,  the  virtual 
work  and  equilibrium  methods  give  the  same  value  for  when 
conventional  finite  difference  expressions  are  used.  The  two 
methods  do  not  give  the  same  result  when  trigonometric  express- 
ions are  used  due  to  the  presence  of  the  two  additional  Fourier 
series  terms  In  the  equilibrium  equation.  Several  major  diff- 
erences were  also  noted  in  the  two  methods.  For  example,  it 
is  more  difficult  to  predict  the  optimum  value  of  A for  the 
equilibrium  approach.  Additionally,  it  was  found  that  an  error 
in  the  estimate  of  A produces  a larger  error  in  the  computed 
value  of  for  the  equilibrium  method  than  for  the  virtual 
work  method.  For  these  reasons,  the  virtual  work  method  is 
recommended  for  general  use  over  the  equilibrium  method.  The 
trigonometric  approach  to  the  finite  difference  calculus  is 

I 

[ recommended  over  the  conventional  approach,  particularly  in 

; those  cases  when  the  shape  of  the  displacement  function  Is 

known  within  rather  broad  tolerances. 
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tan  0 ~ sin  0 ::  0 


cos  0 ~ 1 


ds  = dx 


On  this  basis,  the  equilibrium  equations  are 


EFy  = 0:  qdx  - + (V^+dV^)  = 0 


EM^  =0,  M - PdV  - V^dx  + qdx  ^ - (M+dM)  = 0 


Equations  (A4)  and  (A5)  reduce  to 


V = , M - P iv 
’'f  dx  3x- 


(A6) 

(A7) 


Substitution  of  Eq  (A7)  in  Eq  (A6)  yields 


dx' 


(A8) 


dx 


p p 

Since  M = EI(d  V/dx  ),  the  last  equation  can  be  expressed  as 


d^V  ^ _L  d^V  _ _g_ 
El  7?  "El- 


(A9) 

dx"  dx' 

For  the  case  in  which  the  transverse  loading  Is  zero,  Eq  (A9) 
reduces  to 

(AlO) 

dx^  dx'^ 


116 


Appendix  B 

Development  of  the  Virtual  Work  Equation 

The  shear  and  stress  acting  on  an  element  of  an  arbitrary 
body  with  force,  P,  are  shown  In  Fig.  B.l. 


Fig.  B.l  Element  of  an  arbitrary  body 


The  total  strain  energy  can  be  approximated  by 

U » III  j e dxdydz. 

V 

Substituting  e = a/E  results  In 


(Bl) 


U 


dV., . 


IZ  “''ol 


IB2) 


Since  o * -My/I  and  ^ y^^A,  Eq  (B2)  reduces  to 

U - I ^ dx.  (B3) 

The  Internal  virtual  work  during  buckling  based  on  Eq  {B3)  Is 

6U  - I ^ 6M  dx.  (B4) 
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2 2 

Using  the  expression  M = EI(d  v/dx  ),  Eq  (B4)  can  be  expressed 

6U  = I ^ 6M  dx  (B5) 

The  work  performed  by  the  external  force  on  the  beam  due  to  a 
small  displacement  Is 

Wg  = 1 PA  (B6) 

where  A Is  the  horizontal  displacement  of  the  beam  boundary 
during  buckling.  An  expression  for  the  displacement  of  an 
arbitrary  point,  B,  along  the  beam  can  be  found  from  geometric 
properties  as  shown  In  Fig.  B.2. 


An  expression  for  AB'  Is 

1/2 

AB'  . [dx^  - dx)2]  (B7) 

Using  the  binomial  expansion,  Eq  (B7)  can  be  approximated  as 

AB ' = dx[l  - ^ +•..]  (B8) 

An  expression  for  the  horizontal  displacement  of  B Is 
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AB  - AB'  = ^ dx  (B9) 

The  total  horizontal  displacement  is  found  by  integrating 
Eq  (B9)  along  the  length  of  the  beam: 

A = } I dx  (BIO) 

The  virtual  work  of  the  external  forces  can  be  expressed  as 

6Wg  = P6A  (Bll) 

or 

The  change  in  total  potential  energy,  tv,  can  be  expressed  in 

2 

terms  of  the  variation,  6Tr,  and  the  second  variation,  a tv,  as 
given  by  Taylor's  expansion, 

Air  = Sir  + ^ a^TV  (B1 3) 

It  is  apparent  that  Atv  = 0 since  the  system  is  in  equilibrium. 
For  neutral  equilibrium  (which  corresponds  to  the  onset  of 
buckling),  it  is  necessary  for  a tv  to  also  equal  zero.  Thus 
neutral  equilibrium  corresponds  to  the  case  in  which  Air  = 0, 
The  change  in  total  potential  energy  due  to  a small  virtual 
displacement  is  simply 

Att  = 6U  - 6Wg  (B14) 

Therefore,  the  point  of  neutral  equilibrium  can  be  determined 
from  the  following  condition 
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0 


(B15) 


6U  - 6Wg  = 

Using  this  virtual  work  concept,  Eqs  (B4)  and  (B12)  can  be 
equated  to  give  the  virtual  work  equation: 

El  f dx  - P f ^ ^ dx  = 0 (B16) 
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Appendix  C 

Deri vation  of  Tri gonometri c 
Finite  Difference  Approximati on 

Derivation  of  Second  Derivative 

The  first  three  terms  of  the  Fourier  expansion  can  be 
expressed  as 


V(x)  = 


sin 


tt(x-x  ) it(x-x-) 

— + T3  cos  — ^ 


(Cl) 


Evaluation  of  Eq  (Cl)  at  the  points  x = x^+h  and  x 
suits  in 

= T,  * T2  sin  T3  cos  f 

V.,  = T,  - sin  f * Tj  cos  f 


Xo'h  J-e- 


(C2) 

(C3) 


By  adding  Eqs  (C2)  and  (C3)  and  subtracting  two  times  Eq  (Cl) 


evaluated  at  x = x^.  the  following  expression  is  obtained: 


* ''-1  • 2^3  (cos  I) 


(C4) 


The  second  derivative  of  Eq  (Cl)  with  respect  to  x evaluated  at 

(C5) 


X = Xo  is 


_2  ^(x.-x.) 
Y"  = -T  ~ cos  0 0 
¥0  1 3 ^2  X 


Solving  for  T^  yields 


T.  =-(-)'  V" 
3 7T  0 


(C6) 


If  Eq  (C6)  is  substituted  in  Eq  (C4)  and  the  resulting  express- 
ion solved  for  V",  the  following  equation  is  obtained: 
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or 


2.  . 2 /irh  \ 
4X  sin 


(C7) 


where 


h = 2Xs1  n(iTh/2X)/iT 

Derivation  of  Fourth  Derivative 

The  first  five  terms  of  the  Fourier  serie: 
deriving  the  fourth  derivative  approximation: 


(C8) 

(C9) 


" ^2 

sin  - 

' n ' 

A * ^3 

’’(x-x.) 

s in2 

Evaluation  of  Eq  (CIO)  at  the  points  XQ-2h,  x^-h,  x^,  x^+h , 
and  XQ+2h  results  in  the  following  equations: 


r'2 


X ^ ‘3 


irh  y 

X ■ U 


irh 

X 


V-l  = ^1-^2  cos 


Vo  = T^+T2+T3 


'1 


r'2 

^^2 


Trh 

X 


IIL+  T 

X ^ ' 3 


2A+  T 
'4  ^ 

.2 


Trh 

X 


Trh 

X 


Trh 

X 


1 are  used 

i n 

(x-Xo) 

X 

tt(x-x^) 

X 

(CIO) 

-h,  x„,  x„ 
1 0 0 

+h , 

■5  X 

(Cll) 

cos2 

(C12) 

(C13) 

cosZ  ^ 

(C14) 

;COS4f 

(C15) 

2 12  X 3 X • ig 

Multiplication  of  Eqs  (Cll)  through  (C15)  by  the  appropriate 
factors  and  adding  the  resulting  equations  provides 

''-2  ' ‘>''.1  * - “h  * '‘z  ' '^3  * '^5  - 8T3  cos  f 

^ 8T5  cos2  ^ + 2T3  cos2  ~ + 2T5  cos4  ^ (C16) 
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By  using  trigonometric  substi tuti ons such  as 

cos(20)  = 2(cos0)^-l 
cos(40)  = 8(cos0)^-8(cos0)^+l 


(C17) 

(C18) 


and  combining  terms,  Eq  (C16)  reduces  to 


" "^3  (2cos^-2)^+16T5(sin^^)^  (C19) 


The  fourth  derivative  of  Eq  (CIO)  with  respect  to  x is 


V^''(x)  = (?-^)  [T3+I6T5]. 
A 


(C20) 


Solving  Eq  (C20)  for  and  substituting  the  result  in  Eq  (C19) 
yields  an  expression  for  the  fourth  derivative: 


V'^x)  = (f)'' 


1 


(2COS0-2) 


T (V.2-4V.1+6V, 


4V,+V,)  + 16T.  (1 — K-) 

' 2 ^ (2cos0-2)^ 


(C21) 


where  0 = irh/X,  An  expression  for  Tg  can  be  obtained  by 
applying  Gauss  reduction  to  Eqs  (Cll)  through  (C15).  The 
resulting  expression  is 

_ A,V,.A2V,*A3V„*A2V.,«,V., 


where 


= -2COS0+2 
A2  = 2cos(20)-2 
A3  = -4cos(20)+4cos0 

A^  = 16s1n^0-48cos^0+32cos^0 
+48cos^0-32cos^0 


(C22) 

(C23) 

(C24) 

(C25) 

(C26) 
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In  summary,  the  trigonometric  finite  difference  approximation 
for  the  fourth  derivative  is  given  by  Eq  (C21)  where  Tg  is 
found  by  evaluating  Eq  (C22). 


1 


Half  Station  Trigonometric  Approximation  for  First  Derivative 
The  first  derivative  of  the  Fourier  series  is  shown  by  Eq 


(3-18) 

'''(>■0)  ■ ^2  r 

(3-18) 

and 

(3-19) 

Evaluation  of  Eq  (3-16)  at  XQ+h/2  and  XQ-h/2  results  in 

= T,  + Tj  Sin  T3  cos  (C27) 

'‘-i  ' T,  - Tj  Sin  T3  cos  (C28) 

Subtract  Eq  (C28)  from  Eq  (C27)  to  obtain 

V^l  - = 2T2  sin  (C29) 

If  Eq  (3-19)  is  substituted  in  Eq  (C29)  and  the  terms  rearranged, 
the  following  expression  is  obtained 

V'(x„)  = i (V.,-V  1)  (C30) 

0 j,  +2-2 


where 


I 


h 


2Xsin(|v) 

if 


(C31) 


Appendix  D 
Boundary  Conditions 


The  mathematical  description  of  the  boundary  conditions 
considered  In  this  thesis  can  be  expressed  as  follows: 


1. 

PInned-PInned 

V(0)  » 0 

V(L)  « 0 

V"(0)  = 0 

V(L)  = 0 

2. 

Cl amped-PInned 

V(0)  = 0 

V(L)  = 0 

V'(0)  » 0 

V"(L)  = 0 

3. 

Clamped-Clamped 

V(0)  = 0 

V(L)  » 0 

V'(0)  = 0 

V'(L)  = 0 

4. 

Free-Pinned 

V"(0)  = 0 

V(L)  = 0 

V"  (0}<-PV'  (0)  = 0 

V(L)  = 0 

5. 

Gulded-PInned 

V'(0)  = 0 

V(L)  = 0 

r(o)  = 0 

V''(L)  = 0 

6. 

Clamped-Free 

V(0)  = 0 

V"(L)  = 0 

V'(0)  = 0 

r(L)+TV'(L)  = 

0 

7. 

Cl amped-Gul ded 

V(0)  = 0 

V'(L)  = 0 

V'(0)  = 0 

V'"(L)+7V'(L)  = 

0 

8. 

Gulded-Gulded 

V’(0)  = 0 

V'(L)  = 0 

V"'(0)+'PV' (0)  = 0 

V"'(L)+PV'(L)  = 

0 

9. 

Free-Free 

V"(0)  = 0 

V"(L)  = 0 

V’’'(0)+‘PV’ (0)  = 0 

V"'(L)+PV'(L)  = 

0 

10. 

Free-Guided 

r'(o)+'Fv' (0)  = 0 

V"’(L)+PV'(L)  = 

0 

For  the  virtual 

work  equation,  the 

same  relationship 

) exists 

for  the  virtual  displacements.  For  example,  the  following 
additional  expressions  are  available  for  a pinned-pinned  beam: 
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r ’ 

6V{0)  = 0 

; 6V''(0)  = 0 

[ 6V(L)  = 0 

f. 

I 6V"(L)  = 0 

Finite  difference  expressions  were  substituted  for  the 
derivatives  in  the  above  expressions  to  rel ate  the  external 
and  boundary  grid  points  to  the  internal  grid  points.  In 
I addition,  the  equilibrium  equation  was  applied  at  free  and 

I guided  boundaries  since  there  is  no  end  restraint  for  these 

cases.  Clamped  and  pinned  beams,  on  the  other  hand,  are 
restrained,  and  the  equilibrium  equation  was  not  applied  at 
the  boundary. 

The  third  derivative  boundary  condition  was  not  used  in 
the  virtual  work  approach.  Since  the  boundary  is  free  or 
guided  for  this  boundary  condition,  the  coefficient  of  6V^ 
in  Eq  (3-47)  must  be  zero.  This  additional  equation  provides 
the  necessary  relationship  to  form  the  eigenvalue  problem. 


(Dl) 

(D2) 

(D3) 

(D4) 


Appendix  E 
Summary  of  Resul ts 

Table  II 


Virtual  Work  Approach  with  Five  Grid  Points 


Pinned-Pinned 

1.0 

.75-“ 

-1.3% 

-2.3% 

Cl amped-Pi nned 

.7 

. 5-® 

-4.7% 

-5.7% 

Clamped-Cl amped 

.5 

. 5-® 

-8.1% 

-9.1% 

Free-Pinned 

.75 

.75-“ 

-1.3% 

-2.3% 

Guided-Pinned 

2.0 

1.5  -® 

.44% 

- .57% 

Clamped-Free 

2.0 

1.5  -® 

.44% 

- .57% 

Cl amped-Gui ded 

1.0 

.75-® 

-1.3% 

-2.3% 

Gui ded-Guided 

1.0 

.75-® 

-1.3% 

-2.3% 

Free-Free 

1.0 

.75-® 

-1.3% 

-2.3% 

Free-Guided 

2.0 

1.5  -“ 

.44% 

- .57% 

Table  III 


Equilibrium  Approach  with  Five  Grid  Points 


Boundary 

Condition 

Opt 

X 

Range 
of  X 

Error  for 
X=3.75 

. ■ I.  1 

Conventional 

Error 

Pi nned-Pi nned 

1.0 

1.0 

2.0-~ 

-.9U 

-2.3% 

Cl amped-Pi nned 

.7 

.7 

1.25-“ 

-4.1  % 

-5.7% 

Clamped-Clamped 

.5 

.5 

1 .0-» 

-7.5  % 

-9.1% 

Free-Pinned 

1.0 

1.0 

2.0-» 

-.91% 

-2.3% 

Gui ded-Pi nned 

2.0 

2.0 

3.75-» 

.50% 

- .57% 

Cl amped-Free 

2.0 

2.0 

3.75-» 

.50% 

- .57% 

Clamped-Guided 

1.0 

1.0 

2.0-«o 

-.91% 

-2.3% 

Guided-Gui ded 

1 .0 

1.0 

2.0-“ 

-.91% 

-2.3% 

Free-Free 

1.0 

1.0 

2.0-“ 

-.91% 

-2.3% 

Free-Guided 

i 

2.0 

2.0 

3.75-« 

.50% 

- .57% 
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